
/*
	pbrt source code Copyright(c) 1998-2012 Matt Pharr and Greg Humphreys.

	This file is part of pbrt.

	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions are
	met:

	- Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.

	- Redistributions in binary form must reproduce the above copyright
	notice, this list of conditions and the following disclaimer in the
	documentation and/or other materials provided with the distribution.

	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
	IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
	TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
	PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
	HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
	SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
	LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
	THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
	OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

	*/


// renderers/metropolis.cpp*
#include "stdafx.h"
#include "renderers/metropolis.h"
#include <iostream>

MLTSample::MLTSample(int maxLength) {
	cameraPathSamples.resize(maxLength);
	lightPathSamples.resize(maxLength);
	lightingSamples.resize(maxLength);
}


MLTRenderData::MLTRenderData(MemoryArena &arr, RNG &rand, 
	PathVertex *cp, PathVertex *lp, int cpl, int lpl) : arena(arr), rng(rand), 
	cameraPath(cp), lightPath(lp), cameraPathLength(cpl), lightPathLength(lpl) {
}


MetropolisRenderer::MetropolisRenderer(int perPixelSamples,
	int nboot, int dps, float lsp, bool dds, int mr, int md,
	Camera *c, bool db) {
	camera = c;
	nPixelSamples = perPixelSamples;
	nBootstrap = nboot;
	nDirectPixelSamples = dps;
	maxDepth = md;
	maxConsecutiveRejects = mr;
	bidirectional = db;
	nTasksFinished = 0;
	directLightingIntegrator = dds ?
		new DirectLightingIntegrator(SAMPLE_ALL_UNIFORM, maxDepth) : NULL;

	float largeStepProbability = lsp;
	uint32_t largeStepsPow2 = RoundUpPow2(largeStepProbability * nPixelSamples);
	largeStepsPerPixel = max(1u, largeStepsPow2);
	if (largeStepsPerPixel >= nPixelSamples) largeStepsPerPixel /= 2;
	Assert(largeStepsPerPixel >= 1 && largeStepsPerPixel < nPixelSamples);

	if ((nPixelSamples % largeStepsPerPixel) != 0) {
		int origPixelSamples = nPixelSamples;
		nPixelSamples += largeStepsPerPixel - (nPixelSamples % largeStepsPerPixel);
		Warning("Rounding up to %d Metropolis samples per pixel (from %d)",
			nPixelSamples, origPixelSamples);
	}
}


MetropolisRenderer::~MetropolisRenderer() {
	delete camera;
	delete directLightingIntegrator;
}


MetropolisRenderer *CreateMetropolisRenderer(const ParamSet &params,
	Camera *camera) {
	float largeStepProbability = params.FindOneFloat("largestepprobability", .25f);
	int perPixelSamples = params.FindOneInt("samplesperpixel", 100);
	int nBootstrap = params.FindOneInt("bootstrapsamples", 100000);
	int nDirectPixelSamples = params.FindOneInt("directsamples", 4);
	bool doDirectSeparately = params.FindOneBool("dodirectseparately", true);
	int mr = params.FindOneInt("maxconsecutiverejects", 512);
	int md = params.FindOneInt("maxdepth", 7);
	bool doBidirectional = params.FindOneBool("bidirectional", true);

	if (PbrtOptions.quickRender) {
		perPixelSamples = max(1, perPixelSamples / 4);
		nBootstrap = max(1, nBootstrap / 4);
		nDirectPixelSamples = max(1, nDirectPixelSamples / 4);
	}

	return new MetropolisRenderer(perPixelSamples, nBootstrap,
		nDirectPixelSamples, largeStepProbability, doDirectSeparately,
		mr, md, camera, doBidirectional);
}


void MetropolisRenderer::Render(const Scene *scene) {
	PBRT_MLT_STARTED_RENDERING();
	if (scene->lights.size() > 0) {
		lightDistribution = ComputeLightSamplingCDF(scene);
		vector<PathVertex> cameraPath(maxDepth, PathVertex());
		vector<PathVertex> lightPath(maxDepth, PathVertex());
		MemoryArena arena;
		RNG rng(0);
		MLTRenderData data(arena, rng, &cameraPath[0], &lightPath[0]);
		this->scene = scene;

		if (directLightingIntegrator != NULL) {
			// Compute direct lighting before Metropolis light transport
			LaunchDirectTasks();
			camera->film->WriteImage();
		}

		// Take initial set of samples to compute $b$
		float sumI = 0.f;
		vector<float> bootstrapI;
		float b = Bootstrap(data, &sumI, &bootstrapI);

		// Select initial sample from bootstrap samples
		MLTSample initialSample = GetInitialSample(data, sumI, &bootstrapI[0]);

		// Launch tasks to generate Metropolis samples
		LaunchMLTTasks(data, initialSample, b);
		delete lightDistribution;
	}
	camera->film->WriteImage();
	PBRT_MLT_FINISHED_RENDERING();
}


IntegrationResult MetropolisRenderer::Li(const Scene *scene, 
	const RayDifferential &ray, const Sample *sample, RNG &rng, MemoryArena &arena,
	Intersection *isect, Spectrum *T) const {
	Intersection localIsect;
	if (!isect) isect = &localIsect;
	Spectrum Lo = 0.f;
	if (scene->Intersect(ray, isect))
		Lo = directLightingIntegrator->Li(scene, this, ray, *isect, sample,
		rng, arena).L;
	else {
		// Handle ray that doesn't intersect any geometry
		for (uint32_t i = 0; i < scene->lights.size(); ++i)
			Lo += scene->lights[i]->Le(ray);
	}
	return Lo;
}


Spectrum MetropolisRenderer::Transmittance(const Scene *scene, const RayDifferential &ray,
	const Sample *sample, RNG &rng, MemoryArena &arena) const {
	return 1.f;
}


void MetropolisRenderer::LaunchDirectTasks() {
	PBRT_MLT_STARTED_DIRECTLIGHTING();
	int x0, x1, y0, y1;
	float t0 = camera->shutterOpen, t1 = camera->shutterClose;
	camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);

	if (nDirectPixelSamples > 0) {
		vector<Task *> directTasks;
		LDSampler sampler(x0, x1, y0, y1, nDirectPixelSamples, t0, t1);
		Sample* sample = new Sample(&sampler, directLightingIntegrator, NULL,
			scene);

		int nDirectTasks = max(32 * NumSystemCores(),
			(camera->film->xResolution * camera->film->yResolution) / (16 * 16));
		nDirectTasks = RoundUpPow2(nDirectTasks);

		ProgressReporter directProgress(nDirectTasks, "Direct Lighting");
		for (int i = 0; i < nDirectTasks; ++i)
			directTasks.push_back(new SamplerRendererTask(scene, this, camera,
			directProgress, &sampler, sample, false, i, nDirectTasks));
		std::reverse(directTasks.begin(), directTasks.end());
		EnqueueTasks(directTasks);
		WaitForAllTasks();
		for (uint32_t i = 0; i < directTasks.size(); ++i)
			delete directTasks[i];
		delete sample;
		directProgress.Done();
	}
	PBRT_MLT_FINISHED_DIRECTLIGHTING();
}


void MetropolisRenderer::LaunchMLTTasks(MLTRenderData &data, 
	MLTSample &initialSample, float b) {
	Mutex *filmMutex = Mutex::Create();
	float t0 = camera->shutterOpen, t1 = camera->shutterClose;
	int x0, x1, y0, y1;
	camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);

	uint32_t nTasks = largeStepsPerPixel;
	uint32_t largeStepRate = nPixelSamples / largeStepsPerPixel;
	Info("MLT running %d tasks, large step rate %d", nTasks, largeStepRate);
	Assert(IsPowerOf2(nTasks));

	uint32_t scramble[2] = { data.rng.RandomUInt(), data.rng.RandomUInt() };
	uint32_t pfreq = (x1 - x0) * (y1 - y0);

	ProgressReporter progress(nTasks * largeStepRate, "Metropolis");
	vector<Task *> tasks;
	for (uint32_t i = 0; i < nTasks; ++i) {
		float d[2];
		Sample02(i, scramble, d);
		tasks.push_back(new MLTTask(progress, pfreq, i,
			d[0], d[1], x0, x1, y0, y1, t0, t1, b, initialSample, this, 
			filmMutex));
	}
	EnqueueTasks(tasks);
	WaitForAllTasks();
	for (uint32_t i = 0; i < tasks.size(); ++i)
		delete tasks[i];
	progress.Done();

	Mutex::Destroy(filmMutex);
}


float MetropolisRenderer::Bootstrap(MLTRenderData &data, float *sumI, 
	vector<float> *bootstrapI) {
	PBRT_MLT_STARTED_BOOTSTRAPPING(nBootstrap);
	int x0, x1, y0, y1;
	camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);
	MLTSample sample(maxDepth);

	bootstrapI->reserve(nBootstrap);	
	for (uint32_t i = 0; i < nBootstrap; ++i) {
		// Generate random sample and path radiance for MLT bootstrapping
		float x = Lerp(data.rng.RandomFloat(), x0, x1);
		float y = Lerp(data.rng.RandomFloat(), y0, y1);
		LargeStep(data, &sample, x, y);
		Spectrum L = PathL(data, sample);

		// Compute contribution for random sample for MLT bootstrapping
		float I = calcI(L);
		*sumI += I;
		bootstrapI->push_back(I);
		data.arena.FreeAll();
	}
	float b = *sumI / nBootstrap;
	PBRT_MLT_FINISHED_BOOTSTRAPPING(b);
	Info("MLT computed b = %f", b);
	return b;
}


MLTSample MetropolisRenderer::GetInitialSample(MLTRenderData &data, float sumI,
	float *bootstrapI) {
	float contribOffset = data.rng.RandomFloat() * sumI;
	int x0, x1, y0, y1;
	MLTSample sample(maxDepth);

	camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);
	data.rng.Seed(0);
	sumI = 0.f;

	for (uint32_t i = 0; i < nBootstrap; ++i) {
		float x = Lerp(data.rng.RandomFloat(), x0, x1);
		float y = Lerp(data.rng.RandomFloat(), y0, y1);
		LargeStep(data, &sample, x, y);
		sumI += bootstrapI[i];
		if (sumI > contribOffset)
			break;
	}
	return sample;
}


Spectrum MetropolisRenderer::PathL(MLTRenderData &data, MLTSample &sample) const {
	RayDifferential cameraRay;
	float cameraWt = GenerateCameraRay(sample, &cameraRay);
	data.cameraPathLength = GeneratePath(data, sample, cameraRay, cameraWt, false);

	if (!bidirectional) {
		return Lpath(data, sample);
	}
	else {
		float lightPdf, lightRayPdf;
		
		Ray lightRay;
		Spectrum lightWt;

		lightWt = SampleLightBidir(data, sample, &lightRayPdf, &lightRay);
		if (lightWt.IsBlack() || lightRayPdf == 0.f) {
			return Lpath(data, sample);
		}
		else {
			RayDifferential ray(lightRay);
			data.lightPathLength = GeneratePath(data, sample, ray, lightWt, true);
			return Lbidir(data, sample);
		}
	}
}


float MetropolisRenderer::GenerateCameraRay(const MLTSample &sample,
	RayDifferential *ray) const {
	PBRT_STARTED_GENERATING_CAMERA_RAY((CameraSample *)(&sample.cameraSample));
	float cameraWt = camera->GenerateRayDifferential(sample.cameraSample, ray);
	ray->ScaleDifferentials(1.f / sqrtf(nPixelSamples));
	PBRT_FINISHED_GENERATING_CAMERA_RAY((CameraSample *)(&sample.cameraSample),
		&cameraRay, cameraWt);

	return cameraWt;
}

Spectrum MetropolisRenderer::SampleLightBidir(MLTRenderData &data, 
	MLTSample &sample, float *lightRayPdf, Ray *lightRay) const {
	Normal Nl;
	float lightPdf;
	PBRT_MLT_STARTED_SAMPLE_LIGHT_FOR_BIDIR();
	const float *raySamples = sample.lightRaySamples;

	uint32_t lightNum = lightDistribution->SampleDiscrete(
		sample.lightNumSample, &lightPdf);
	const Light *light = scene->lights[lightNum];
	LightSample lrs(raySamples[0], raySamples[1], raySamples[2]);
	Spectrum lightWt = light->Sample_L(scene, lrs, raySamples[3], raySamples[4],
		sample.cameraSample.time, lightRay, &Nl, lightRayPdf);
	lightWt *= AbsDot(Normalize(Nl), lightRay->d) / (*lightRayPdf * (lightPdf));

	PBRT_MLT_FINISHED_SAMPLE_LIGHT_FOR_BIDIR();

	return lightWt;
}


Spectrum MetropolisRenderer::Lpath(MLTRenderData &data, MLTSample &sample) const {
	PBRT_MLT_STARTED_LPATH();
	Spectrum L(0.f);
	bool prevSpecular = true, allSpecular = true;
	for (int i = 0; i < data.cameraPathLength; ++i) {
		// Initialize basic variables for camera path vertex
		const PathVertex &vc = data.cameraPath[i];
		const Point &pc = vc.bsdf->dgShading.p;
		const Normal &nc = vc.bsdf->dgShading.nn;

		// Add emitted light from vertex if appropriate
		if (prevSpecular && (directLightingIntegrator == NULL || !allSpecular))
			L += vc.alpha * vc.isect.Le(vc.wPrev);

		// Add direct illumination for Metropolis path vertex
		if (directLightingIntegrator == NULL || !allSpecular) {
			L += DirectLightingVertex(data, sample, i, vc, pc, nc);
		}

		prevSpecular = vc.specularBounce;
		allSpecular &= prevSpecular;
	}

	// Add contribution of escaped ray, if any
	L += EscapedRayL(data, prevSpecular, allSpecular);
	PBRT_MLT_FINISHED_LPATH();
	return L;
}


Spectrum MetropolisRenderer::DirectLightingVertex(MLTRenderData &data, 
	MLTSample &sample, int pathIdx, const PathVertex &vc, const Point &pc, 
	const Normal &nc) const {
	Spectrum Ld(0.f);

	// Choose light and call _EstimateDirect()_ for Metropolis vertex
	const LightingSample &ls = sample.lightingSamples[pathIdx];
	float lightPdf;
	uint32_t lightNum = lightDistribution->SampleDiscrete(ls.lightNum,
		&lightPdf);
	const Light *light = scene->lights[lightNum];

	PBRT_MLT_STARTED_ESTIMATE_DIRECT();
	Ld = vc.alpha *
		EstimateDirect(scene, this, data.arena, light, pc, nc, vc.wPrev,
		vc.isect.rayEpsilon, sample.cameraSample.time, vc.bsdf, data.rng,
		ls.lightSample, ls.bsdfSample,
		BxDFType(BSDF_ALL & ~BSDF_SPECULAR)) / lightPdf;
	PBRT_MLT_FINISHED_ESTIMATE_DIRECT();
	return Ld;
}


Spectrum MetropolisRenderer::EscapedRayL(MLTRenderData &data, bool prevSpecular,
	bool allSpecular) const {
	Spectrum L(0.f);
	if (!data.escapedAlpha.IsBlack() && prevSpecular &&
		(directLightingIntegrator == NULL || !allSpecular))
	for (uint32_t i = 0; i < scene->lights.size(); ++i)
		L += (data.escapedAlpha) * scene->lights[i]->Le(data.escapedRay);
	return L;
}


Spectrum MetropolisRenderer::Lbidir(MLTRenderData &data, MLTSample &sample) const {
	PBRT_MLT_STARTED_LBIDIR();
	Spectrum L(0.f);
	bool prevSpecular = true, allSpecular = true;

	// Compute number of specular vertices for each path length
	int nVerts = data.cameraPathLength + data.lightPathLength + 2;
	int *nSpecularVertices = new int[nVerts];
	GetNumSpecularVertices(data, nSpecularVertices);

	for (int i = 0; i < data.cameraPathLength; ++i) {
		// Initialize basic variables for camera path vertex
		const PathVertex &vc = data.cameraPath[i];
		const Point &pc = vc.bsdf->dgShading.p;
		const Normal &nc = vc.bsdf->dgShading.nn;

		// Add emitted light from vertex if appropriate
		if (prevSpecular && (directLightingIntegrator == NULL || !allSpecular))
			L += vc.alpha * vc.isect.Le(vc.wPrev);

		// Compute direct illumination for Metropolis path vertex
		Spectrum Ld(0.f);
		if (directLightingIntegrator == NULL || !allSpecular) {
			Ld = DirectLightingVertex(data, sample, i, vc, pc, nc);
		}

		prevSpecular = vc.specularBounce;
		allSpecular &= prevSpecular;
		L += Ld / (i + 1 - nSpecularVertices[i + 1]);

		if (!vc.specularBounce) {
			// Loop over light path vertices and connect to camera vertex
			for (int j = 0; j < data.lightPathLength; ++j) {
				L += LightPathL(data, sample, pc, vc, nc, nSpecularVertices, i, j);
			}
		}
	}
	// Add contribution of escaped ray, if any
	L += EscapedRayL(data, prevSpecular, allSpecular);
	PBRT_MLT_FINISHED_LBIDIR();

	delete[] nSpecularVertices;
	return L;
}


void MetropolisRenderer::GetNumSpecularVertices(MLTRenderData &data,
	int *nSpecularVertices) const {
	int nVerts = data.cameraPathLength + data.lightPathLength + 2;
	memset(nSpecularVertices, 0, nVerts * sizeof(int));
	for (int i = 0; i < data.cameraPathLength; ++i)
	for (int j = 0; j < data.lightPathLength; ++j)
	if (data.cameraPath[i].specularBounce || data.lightPath[j].specularBounce)
		++nSpecularVertices[i + j + 2];
}


Spectrum MetropolisRenderer::LightPathL(MLTRenderData &data, MLTSample &sample, 
	const Point &pc, const PathVertex &vc, const Normal &nc, 
	int *nSpecularVertices, int cameraPathIdx, int lightPathIdx) const {
	const PathVertex &vl = data.lightPath[lightPathIdx];
	const Point &pl = vl.bsdf->dgShading.p;
	const Normal &nl = vl.bsdf->dgShading.nn;
	Spectrum L(0.f);

	if (!vl.specularBounce) {
		// Compute contribution between camera and light vertices
		Vector w = Normalize(pl - pc);
		Spectrum fc = vc.bsdf->f(vc.wPrev, w) * (1 + vc.nSpecularComponents);
		Spectrum fl = vl.bsdf->f(-w, vl.wPrev) * (1 + vl.nSpecularComponents);
		Ray r(pc, pl - pc, 1e-3f, .999f, sample.cameraSample.time);

		if (!fc.IsBlack() && !fl.IsBlack() && !scene->IntersectP(r)) {
			// Compute weight for bidirectional path, _pathWt_
			int pathIdx = cameraPathIdx + lightPathIdx + 2;
			float pathWt = 1.f / (pathIdx - nSpecularVertices[pathIdx]);
			float G = AbsDot(nc, w) * AbsDot(nl, w) / DistanceSquared(pl, pc);

			L += (vc.alpha * fc * G * fl * vl.alpha) * pathWt;
		}
	}

	return L;
}


void MetropolisRenderer::LargeStep(MLTRenderData &data, MLTSample *sample, 
	float x, float y) const {
	float t0 = camera->shutterOpen, t1 = camera->shutterClose;
	RNG &rng = data.rng;

	// Do large step mutation of _cameraSample_
	sample->cameraSample.imageX = x;
	sample->cameraSample.imageY = y;
	sample->cameraSample.time = Lerp(rng.RandomFloat(), t0, t1);
	sample->cameraSample.lensU = rng.RandomFloat();
	sample->cameraSample.lensV = rng.RandomFloat();
	for (int i = 0; i < maxDepth; ++i) {
		// Apply large step to $i$th camera _PathSample_
		PathSample &cps = sample->cameraPathSamples[i];
		cps.bsdfSample.uComponent = rng.RandomFloat();
		cps.bsdfSample.uDir[0] = rng.RandomFloat();
		cps.bsdfSample.uDir[1] = rng.RandomFloat();
		cps.rrSample = rng.RandomFloat();

		// Apply large step to $i$th _LightingSample_
		LightingSample &ls = sample->lightingSamples[i];
		ls.bsdfSample.uComponent = rng.RandomFloat();
		ls.bsdfSample.uDir[0] = rng.RandomFloat();
		ls.bsdfSample.uDir[1] = rng.RandomFloat();
		ls.lightNum = rng.RandomFloat();
		ls.lightSample.uComponent = rng.RandomFloat();
		ls.lightSample.uPos[0] = rng.RandomFloat();
		ls.lightSample.uPos[1] = rng.RandomFloat();
	}
	if (bidirectional) {
		// Apply large step to bidirectional light samples
		sample->lightNumSample = rng.RandomFloat();
		for (int i = 0; i < 5; ++i)
			sample->lightRaySamples[i] = rng.RandomFloat();
		for (int i = 0; i < maxDepth; ++i) {
			// Apply large step to $i$th light _PathSample_
			PathSample &lps = sample->lightPathSamples[i];
			lps.bsdfSample.uComponent = rng.RandomFloat();
			lps.bsdfSample.uDir[0] = rng.RandomFloat();
			lps.bsdfSample.uDir[1] = rng.RandomFloat();
			lps.rrSample = rng.RandomFloat();
		}
	}
}


void MetropolisRenderer::mutate(RNG &rng, float *v, float min, float max) const {
	if (min == max) { *v = min; return; }
	Assert(min < max);
	float a = 1.f / 1024.f, b = 1.f / 64.f;
	static const float logRatio = -logf(b / a);
	float delta = (max - min) * b * expf(logRatio * rng.RandomFloat());
	if (rng.RandomFloat() < 0.5f) {
		*v += delta;
		if (*v >= max) *v = min + (*v - max);
	}
	else {
		*v -= delta;
		if (*v < min) *v = max - (min - *v);
	}
	if (*v < min || *v >= max) *v = min;
}


void MetropolisRenderer::SmallStep(MLTRenderData &data, MLTSample *sample) const {
	float t0 = camera->shutterOpen, t1 = camera->shutterClose;
	RNG &rng = data.rng;
	int x0, x1, y0, y1;
	camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);

	mutate(rng, &sample->cameraSample.imageX, x0, x1);
	mutate(rng, &sample->cameraSample.imageY, y0, y1);
	mutate(rng, &sample->cameraSample.time, t0, t1);
	mutate(rng, &sample->cameraSample.lensU);
	mutate(rng, &sample->cameraSample.lensV);
	// Apply small step mutation to camera, lighting, and light samples
	for (int i = 0; i < maxDepth; ++i) {
		// Apply small step to $i$th camera _PathSample_
		PathSample &eps = sample->cameraPathSamples[i];
		mutate(rng, &eps.bsdfSample.uComponent);
		mutate(rng, &eps.bsdfSample.uDir[0]);
		mutate(rng, &eps.bsdfSample.uDir[1]);
		mutate(rng, &eps.rrSample);

		// Apply small step to $i$th _LightingSample_
		LightingSample &ls = sample->lightingSamples[i];
		mutate(rng, &ls.bsdfSample.uComponent);
		mutate(rng, &ls.bsdfSample.uDir[0]);
		mutate(rng, &ls.bsdfSample.uDir[1]);
		mutate(rng, &ls.lightNum);
		mutate(rng, &ls.lightSample.uComponent);
		mutate(rng, &ls.lightSample.uPos[0]);
		mutate(rng, &ls.lightSample.uPos[1]);
	}

	if (bidirectional) {
		mutate(rng, &sample->lightNumSample);
		for (int i = 0; i < 5; ++i)
			mutate(rng, &sample->lightRaySamples[i]);
		for (int i = 0; i < maxDepth; ++i) {
			// Apply small step to $i$th light _PathSample_
			PathSample &lps = sample->lightPathSamples[i];
			mutate(rng, &lps.bsdfSample.uComponent);
			mutate(rng, &lps.bsdfSample.uDir[0]);
			mutate(rng, &lps.bsdfSample.uDir[1]);
			mutate(rng, &lps.rrSample);
		}
	}
}

// Metropolis Method Definitions
uint32_t MetropolisRenderer::GeneratePath(MLTRenderData &data, MLTSample &sample,
	const RayDifferential &r, const Spectrum &a, bool isLightPath) const {
	PBRT_MLT_STARTED_GENERATE_PATH();
	RayDifferential ray(r);
	Spectrum alpha(a);

	vector<PathSample> *samples;
	PathVertex *path = isLightPath ? data.lightPath : data.cameraPath;
	samples = isLightPath ? &sample.lightPathSamples : &sample.cameraPathSamples;
	if (!isLightPath) data.escapedAlpha = 0.f;

	uint32_t length = 0;
	for (; length < samples->size(); ++length) {
		// Try to generate next vertex of ray path
		PathVertex *v = &path[length];
		if (!GeneratePathVertex(data, ray, alpha, isLightPath, v)) break;

		// Sample direction for outgoing Metropolis path direction
		float pdf;
		BxDFType flags;
		Spectrum f = v->bsdf->Sample_f(-ray.d, &v->wNext, 
			(*samples)[length].bsdfSample, &pdf, BSDF_ALL, &flags);
		v->specularBounce = (flags & BSDF_SPECULAR) != 0;
		v->nSpecularComponents = v->bsdf->NumComponents(BxDFType(BSDF_SPECULAR |
			BSDF_REFLECTION | BSDF_TRANSMISSION));
		if (f.IsBlack() || pdf == 0.f || 
			!RussianRoulette(*v, f, pdf, 
			(*samples)[length].rrSample, &ray, &alpha))
		{
			PBRT_MLT_FINISHED_GENERATE_PATH();
			return length + 1;
		}
	}
	PBRT_MLT_FINISHED_GENERATE_PATH();
	return length;
}


bool MetropolisRenderer::GeneratePathVertex(MLTRenderData &data, 
	RayDifferential &ray, Spectrum &alpha, bool isLightPath, PathVertex *v) const {
	if (!scene->Intersect(ray, &v->isect)) {
		// Handle ray that leaves the scene during path generation
		if (!isLightPath) {	// TODO double check negation
			data.escapedAlpha = alpha;
			data.escapedRay = ray;
		}
		return false;
	}

	// Record information for current path vertex
	v->alpha = alpha;
	BSDF *bsdf = v->isect.GetBSDF(ray, data.arena);
	v->bsdf = bsdf;
	v->wPrev = -ray.d;

	return true;
}


bool MetropolisRenderer::RussianRoulette(PathVertex &v, Spectrum &f, float pdf, 
	float rrSample, RayDifferential *ray, Spectrum *alpha) const {
	// Terminate path with RR or prepare for finding next vertex
	const Point &p = v.bsdf->dgShading.p;
	const Normal &n = v.bsdf->dgShading.nn;
	Spectrum pathScale = f * AbsDot(v.wNext, n) / pdf;
	float rrSurviveProb = min(1.f, pathScale.y());

	if (rrSample > rrSurviveProb) return false;

	*alpha *= pathScale / rrSurviveProb;
	//alpha *= renderer->Transmittance(scene, ray, NULL, rng, arena);
	*ray = RayDifferential(p, v.wNext, *ray, v.isect.rayEpsilon);

	return true;
}


MLTTask::MLTTask(ProgressReporter &prog, uint32_t pfreq, uint32_t tn,
	float ddx, float ddy, int xx0, int xx1, int yy0, int yy1, float tt0, 
	float tt1, float bb, const MLTSample &is, MetropolisRenderer *ren, Mutex *fm)
	: progress(prog), initialSample(is) {
	progressUpdateFrequency = pfreq;
	taskNum = tn;
	dx = ddx;
	dy = ddy;
	x0 = xx0;
	x1 = xx1;
	y0 = yy0;
	y1 = yy1;
	t0 = tt0;
	t1 = tt1;
	currentPixelSample = 0;
	b = bb;
	renderer = ren;
	filmMutex = fm;
}

void MLTTask::Run() {
	PBRT_MLT_STARTED_MLT_TASK(this);
	// Declare variables for storing and computing MLT samples
	vector<PathVertex> cameraPath(renderer->maxDepth, PathVertex());
	vector<PathVertex> lightPath(renderer->maxDepth, PathVertex());
	MemoryArena arena;
	RNG rng(taskNum);
	MLTRenderData data(arena, rng, &cameraPath[0], &lightPath[0]);

	// Declare basic _MLTTask_ variables and prepare for sampling
	vector<MLTSample> samples(2, MLTSample(renderer->maxDepth));
	Spectrum L[2];
	float I[2];
	Init(data, &samples[0], L, I);

	for (uint64_t s = 0; s < nTaskSamples; ++s) {
		// Compute proposed mutation to current sample
		samples[proposed] = samples[current];
		Mutate(data, s, &samples[proposed]);

		// Compute contribution of proposed sample
		L[proposed] = renderer->PathL(data, samples[proposed]);
		I[proposed] = renderer->calcI(L[proposed]);
		data.arena.FreeAll();

		// Compute acceptance probability for proposed sample
		float acceptProb = min(1.f, I[proposed] / I[current]);

		// Splat current and proposed samples to _Film_
		SplatSample(acceptProb, &samples[0], L, I);

		// Randomly accept proposed path mutation (or not)
		VetMutation(data, acceptProb);
	}
	Assert(pixelNumOffset == nPixels);

	// Update display for recently computed Metropolis samples
	UpdateDisplay();

	PBRT_MLT_FINISHED_MLT_TASK(this);
}

void MLTTask::Init(MLTRenderData &data, MLTSample *samples, Spectrum *L, float *I) 
{
	PBRT_MLT_STARTED_TASK_INIT();
	nPixels = (x1 - x0) * (y1 - y0);
	nPixelSamples = renderer->nPixelSamples;
	largeStepRate = nPixelSamples / renderer->largeStepsPerPixel;
	Assert(largeStepRate > 1);
	nTaskSamples = uint64_t(nPixels) * uint64_t(largeStepRate);
	consecutiveRejects = 0;
	progressCounter = progressUpdateFrequency;

	current = 0;
	proposed = 1;

	// Compute _L[current]_ for initial sample
	samples[current] = initialSample;
	L[current] = renderer->PathL(data, samples[current]);
	I[current] = renderer->calcI(L[current]);
	data.arena.FreeAll();

	// Compute randomly permuted table of pixel indices for large steps
	pixelNumOffset = 0;
	largeStepPixelNum.reserve(nPixels);
	for (uint32_t i = 0; i < nPixels; ++i) largeStepPixelNum.push_back(i);
	Shuffle(&largeStepPixelNum[0], nPixels, 1, data.rng);
	PBRT_MLT_FINISHED_TASK_INIT();
}

void MLTTask::Mutate(MLTRenderData &data, uint64_t sampleIdx, MLTSample *sample) {
	PBRT_MLT_STARTED_MUTATION();
	bool largeStep = ((sampleIdx % largeStepRate) == 0);
	if (largeStep) {
		int x = x0 + largeStepPixelNum[pixelNumOffset] % (x1 - x0);
		int y = y0 + largeStepPixelNum[pixelNumOffset] / (x1 - x0);
		renderer->LargeStep(data, sample, x + dx, y + dy);
		++pixelNumOffset;
	}
	else
		renderer->SmallStep(data, sample);
	PBRT_MLT_FINISHED_MUTATION();
}

void MLTTask::SplatSample(float acceptProb, MLTSample *samples, Spectrum *L, 
	float *I) {
	PBRT_MLT_STARTED_SAMPLE_SPLAT();
	if (I[current] > 0.f) {
		if (!isinf(1.f / I[current])) {
			Spectrum contrib = (b / nPixelSamples) * L[current] / I[current];
			renderer->camera->film->Splat(samples[current].cameraSample, 
				Spectrum((1.f - acceptProb) * contrib));
		}
	}
	if (I[proposed] > 0.f) {
		if (!isinf(1.f / I[proposed])) {
			Spectrum contrib = (b / nPixelSamples) * L[proposed] / I[proposed];
			renderer->camera->film->Splat(samples[proposed].cameraSample, 
				Spectrum(acceptProb * contrib));
		}
	}
	PBRT_MLT_FINISHED_SAMPLE_SPLAT();
}

void MLTTask::VetMutation(MLTRenderData &data, float acceptProb) {
	if (consecutiveRejects >= renderer->maxConsecutiveRejects ||
		data.rng.RandomFloat() < acceptProb) {
		PBRT_MLT_ACCEPTED_MUTATION(
			acceptProb, &samples[current], &samples[proposed]);
		current ^= 1;
		proposed ^= 1;
		consecutiveRejects = 0;
	}
	else
	{
		PBRT_MLT_REJECTED_MUTATION(
			acceptProb, &samples[current], &samples[proposed]);
		++consecutiveRejects;
	}
	if (--progressCounter == 0) {
		progress.Update();
		progressCounter = progressUpdateFrequency;
	}
}


void MLTTask::UpdateDisplay() {
	PBRT_MLT_STARTED_DISPLAY_UPDATE();
	int ntf = AtomicAdd(&renderer->nTasksFinished, 1);
	int64_t totalSamples = int64_t(nPixels) * int64_t(nPixelSamples);
	float splatScale = float(double(totalSamples) / double(ntf * nTaskSamples));
	renderer->camera->film->UpdateDisplay(x0, y0, x1, y1, splatScale);
	if ((taskNum % 8) == 0) {
		MutexLock lock(*filmMutex);
		renderer->camera->film->WriteImage(splatScale);
	}
	PBRT_MLT_FINISHED_DISPLAY_UPDATE();
}